Before using R to illustrate basic programming concepts and data analysis tools, we will get familiar with the RStudio layout.
RStudio has four primary panels that will help you interact with your data. We will use the default layout of these panels.
Let’s use these panels to create and interact with data.
Console:
2 + 2 into the console panel and hit ENTERsum = 2 + 2 into the console panel and hit ENTERSource:
sum = 2 + 3 into the source panel and hit cntrl+ENTEREnvironment:
sum is stored in our environmentrm(sum) to clear the object from the environmentViewer:
Now that we understand the layout, we are ready to review the concepts covered in Module 2 Week 2.2. These concepts will help us understand what is happening when we create and manipulate data.
“Object” is a generic term for anything that R stores in the environment. This can include anything from an individual number or word, to lists of values, to entire datasets.
Importantly, objects belong to different “classes” depending on the type of values that they store.
"hello world" and "welcome to R".TRUE or FALSE# Create a numeric object
my_number = 5.6
# Check the class
class(my_number)[1] "numeric"# Create a character object
my_character = "welcome to R"
# Check the class
class(my_character)[1] "character"# Create a logical object
my_logical = TRUE
# Check the class
class(my_logical)[1] "logical"R can perform operations on objects.
# Create a numeric object
my_number = 5.6
# Check the class
class(my_number)[1] "numeric"# Perform a calculation
my_number = my_number + 5The class of an object determines the type of operations you can perform on it. Some operations can only be run on numeric objects (numbers).
# Create a character object
my_number = "5.6"
# Check the class
class(my_number)
# Perform a calculation
my_number + 5
round(my_number)R contains functions that can convert some objects to different factors.
# Convert character to numeric
my_number = as.numeric("5")
class(my_number)[1] "numeric"my_number <- 5
# But R is only so smart
my_number = as.numeric("five")
print(my_number)[1] NAThe most simple objects are single values, but most data analysis involves more complicated data structures.
Lists are a type of data structure that store multiple values together. Lists are created using c() and allow you to perform operations on a series of values.
# Create a numeric list (also called a "vector")
numeric_vector = c(6, 11, 13, 31)
# Print the vector
print(numeric_vector)[1] 6 11 13 31# Check the class
class(numeric_vector)[1] "numeric"# Calculate the mean
mean(numeric_vector)[1] 15.25An important part of working with more complex data structures is called “indexing.” Indexing allows you to extract specific values from a data structure.
# Extract the 2nd element from the list
numeric_vector[2][1] 11# Extract elements 2-4
numeric_vector[2:4][1] 11 13 31# Extract elements 1-2
numeric_vector[c(TRUE, TRUE, FALSE, FALSE)][1] 6 11Data frames are the most common type of data structure used in research. Data frames combine multiple lists of values into a single object.
# Create a dataframe
my_data = data.frame(
x1 = rnorm(100, mean = 1, sd = 1),
x2 = rnorm(100, mean = 1, sd = 1)
)
class(my_data)[1] "data.frame"Anything that comes in a spreadsheet (for example, an excel file) can be loaded into an R environment as a dataframe. R works most easily when spreadsheets are saved as a .csv file.
# Use `read.csv()` to load data from a website
dat = read.csv("https://raw.githubusercontent.com/jrspringman/psci3200-globaldev/main/workshops/aau_survey/clean_endline_did.csv")
# Use `read.csv()` to load data from your computer's Downloads folder
# dat = read.csv("/home/jeremy/Downloads/clean_endline_did.csv")In most data frames, rows correspond to observations and the columns correspond to variables that describe the observations. Here, we are looking at survey data from an RCT involving university students in Addis Ababa. Each row correspondents to a different survey respondent, and each column represents their answers to a different question from the survey.
Packages are an extremely important part of data analysis with R.
install.packages() to install a package and library() to load a packageIn the next section, we’ll use the package dplyr to perform some data cleaning. dplyr is part of a universe of packages called tidyverse. Since this is one of the most important packages in the R ecosystem, let’s install and load it.
# install.packages("tidyverse")
library(tidyverse)When you are searching online or asking ChatGPT how to perform a specific task in R, it often helps to specify that you are looking for a solution in dplyr.
In the real-world, data never comes ready to be analyzed. Data cleaning is the process of manipulating data so that it can be analyzed. This is usually the most difficult and time-consuming part of any data analysis project. Let’s walk through some examples.
Question: In your work, what data cleaning have you had to do? How have you cleaned your data?
Imagine we want to analyze the relationship between whether a respondent moved to to Addis Ababa to attend university and their level of political participation. However, there are two problems:
How can we create a variable measuring whether the respondent moved to Addis Ababa? We have a multiple-choice question asking students about what region they come from.
Let’s start by investigating this variable.
# The name of the variable in our dataframe is `q8_baseline`
table(dat$q8_baseline)
Addis Ababa
283
Afar Region
1
Amhara Region
174
Benishangul-Gumuz Region
5
Dire Dawa (city)
5
Gambela Region
1
Harari Region
4
Oromia Region
205
Other
12
Prefer not to say
19
Sidama Region
35
Somali Region
5
South West Ethiopia Peoples Region
6
Southern Nations, Nationalities, and Peoples Region
65
Tigray Region
5 dat = dat %>% # this is called a "pipe"
# give our variable a better name
rename(home_region = q8_baseline)
dat = dat %>%
# drop respondents who report "Prefer not to say"
filter(!home_region == "Prefer not to say")
dat = dat %>%
# clean home region variable using `mutate()`
mutate(
# Shorten a long name to an abbreviation
home_region = ifelse(home_region == "Southern Nations, Nationalities, and Peoples Region", "SNNPR", home_region),
# remove the word "Region" from every observation of this column
home_region = str_remove(home_region, " Region"),
home_region = str_remove(home_region, " \\(city\\)")
)
# Check if it worked
table(dat$home_region)
Addis Ababa Afar
283 1
Amhara Benishangul-Gumuz
174 5
Dire Dawa Gambela
5 1
Harari Oromia
4 205
Other Sidama
12 35
SNNPR Somali
65 5
South West Ethiopia Peoples Tigray
6 5 # Chain these all together for more concise code
dat = read_csv("https://raw.githubusercontent.com/jrspringman/psci3200-globaldev/main/workshops/aau_survey/clean_endline_did.csv" ) %>%
# give our variable a better name
rename(home_region = q8_baseline) %>%
# drop respondents who report "Prefer not to say"
filter(!home_region == "Prefer not to say") %>%
# clean home region variable
mutate(
# Shorten a long name to an abbreviation
home_region = case_when(
home_region == "Southern Nations, Nationalities, and Peoples Region" ~ "SNNPR",
home_region == "South West Ethiopia Peoples Region" ~ "SWEPR",
TRUE ~ home_region
),
# remove the word "Region" from every observation of this column
home_region = str_remove(home_region, " Region| \\(city\\)")
)
# Check if it worked
table(dat$home_region)
Addis Ababa Afar Amhara Benishangul-Gumuz
283 1 174 5
Dire Dawa Gambela Harari Oromia
5 1 4 205
Other Sidama SNNPR Somali
12 35 65 5
SWEPR Tigray
6 5 Now that we’ve cleaned-up the names, we want to create a variable that tells us whether or not each respondent is originally from Addis Ababa. This will let us measure whether or not they moved to Addis in order to attend college.
# Creating a measure of whether a respondent moved to Addis Ababa in order to at independent variable
dat = dat %>%
mutate(
moved = case_when(
home_region == "Addis Ababa" ~ 0,
TRUE ~ 1
)
)
table(dat$moved, dat$home_region)
Addis Ababa Afar Amhara Benishangul-Gumuz Dire Dawa Gambela Harari Oromia
0 283 0 0 0 0 0 0 0
1 0 1 174 5 5 1 4 205
Other Sidama SNNPR Somali SWEPR Tigray
0 0 0 0 0 0 0
1 12 35 65 5 6 5Now, we need to create our second variable measuring levels of political participation. Remember, the challenge is that we have multiple measures of participation. Let’s start with two measures:
q13_4_1q13_5_1# Check out the distribution of our variables
dat %>% select(q13_4_1, q13_5_1) %>% head(5)# A tibble: 5 × 2
q13_4_1 q13_5_1
<dbl> <dbl>
1 0 0
2 3 0
3 0 0
4 0 0
5 0 2mean(dat$q13_4_1)[1] NAhist(dat$q13_5_1)
# Create a single measuree
dat = dat %>%
mutate(
add_participation = q13_4_1 + q13_5_1
)
hist(dat$add_participation)
# how do investigate NA values?One thing we need to be careful with are NA values. We need to think carefully about why NA values are in our data and how to handle them appropriately.
Question: Thinking about the data that you have worked with, what are the most common sources of NA values?
## Find our two participation measures
add_ecols = grep("q13_4_1$|q13_5_1$", names(dat), value = T)
dat = dat %>%
mutate(add_participation = rowSums(across(add_ecols) ) )
#mutate(add_participation_end = rowSums(across(add_ecols), na.rm = T) )
# dat = dat %>%
# mutate(home_region = na_if(home_region, "Prefer not to say")) %>%Often, data analysis projects will require you to use variables from more than one dataset. This will require you to combine separate datasets into a single dataset in R, which is called merging. A merge uses an identifier (administrative unit names, respondent IDs, etc.) that is present in both datasets to combine them into one.
Question: Thinking about the data that you have worked with, have you ever had to merge multiple datasets? What was the data? How did you do it?
## Here we are selecting the unit-identifier (response_id) and the two variables we have created and saving them to a separate dataframe
new_vars = dat %>%
select(response_id, moved, add_participation)
## After you do some cleaning or create new variables, you may want to save that dataframe for future use
# write.csv(new_dat, here::here("workshops/aau_survey/new_vars.csv"))
## Drop the new variables from our dataframe so that we can merge them back in
dat = dat %>%
select( -moved, -add_participation)
## Check to make sure the merge will work as expected
# nrow(dat)
# nrow(new_vars)
# dat$response_id[! dat$response_id %in% new_vars$response_id]
## Use dplyr's `full_join` function to merge datasets
dat = dat %>%
full_join(., new_vars)One of the major benefits of using R is the ability to create beautiful tables and figures that can be incorporated into deliverables and other professional documents. Importantly, you can create a script that allows you to automatically recreate the same visualizations when you want to make slight changes or incorporate new data.
Question: What methods do you currently use to create data visualizations for deliverables and other research products?
#install.packages("gt")
library(gt)
desciptives = dat %>%
group_by(moved) %>%
summarise(observations = n(),
share = observations / nrow(dat),
mean = mean(add_participation, na.rm = T)) %>%
mutate(moved = case_when(moved == 0 ~ "No",
moved == 1 ~ "Yes"),
share = share *100)
gt(desciptives) %>%
tab_header(
title = "Moving to University and Political Participation",
) %>%
fmt_number(
columns = 3:4,
decimals = 2,
use_seps = FALSE
) %>%
cols_label(
moved = md("Moved to<br>Addis"),
observations = md("Respondents<br>(count)"),
share = md("Respondents<br>(%)"),
mean = md("Participation<br>(mean score)")
) %>%
gtsave(filename = "tab_1.png")There are many ways to create figures in R, but ggplot is the most popular.
ggplot(dat, aes(add_participation)) +
geom_histogram(binwidth = 1)
dat %>%
filter(add_participation < 15) %>%
ggplot(aes(add_participation)) +
geom_histogram(binwidth = 1)
dat %>%
filter(add_participation < 15) %>%
ggplot(aes(x = add_participation, fill = q3_baseline)) +
geom_histogram(binwidth = 1, alpha = 0.5, position = "identity")
dat %>%
filter(add_participation < 15) %>%
mutate(moved = fct_rev(as.factor(moved))) %>%
ggplot(aes(x = add_participation, fill = moved)) +
geom_histogram(aes(y = (..count..) / sum(..count..) * 100),
binwidth = 1, alpha = 0.5, position = "identity") +
scale_y_continuous(labels = scales::percent_format(scale = 1)) +
labs(y = "Percentage") +
theme_bw() +
theme(legend.title = element_blank(), legend.position = c(.8, .8))
You can do tremendous amounts of customization with ggplot to create extremely informative and professional plots.
## demographics
demos = dat %>%
drop_na(q3_baseline) %>%
select(`Respondent Gender` = q3_baseline, `Work as a student?` = q4_baseline,
`Rural or urban?` = q5_baseline, `Financial support from family?` = q6_baseline,
`Home region` = home_region,
`Student Year` = class_year) %>%
pivot_longer(everything()) %>%
group_by(name, value) %>%
tally() %>%
mutate(pct = n/sum(n)) %>%
top_n(n = 5, wt = pct)
demos$name = factor(demos$name, levels = c('Home region', 'Rural or urban?', 'Student Year', 'Respondent Gender', 'Work as a student?', 'Financial support from family?'))
demos <- demos %>%
group_by(name) %>%
mutate(total_n = sum(n)) %>%
ungroup()
ggplot(demos , aes(y = value, x = pct)) +
geom_col(fill = "grey") +
facet_wrap( ~name, scales = "free") +
scale_y_discrete(labels = scales::label_wrap(30)) +
hrbrthemes::scale_x_percent(limits = c(0, 1)) +
labs(x = "Percent of respondents to the survey", y = NULL,
title = "Demographic characteristics of baseline respondents",
subtitle = glue::glue("Number of respondents = {scales::comma(nrow(df))}"),
caption = "Note: home region only displays top five categories by size.") +
geom_text(
aes(
label = glue::glue("n = {total_n}")
),
x = 0.92, y = -Inf,
vjust = -1, hjust = 1,
inherit.aes = FALSE
)
#ggsave("/home/jeremy/Downloads/demographics.png")Now let’s use simple linear regression to estimate the relationship between these two variables that we created.
Question: Have you used regression in your research? If so, what software did you use to
est = lm(add_participation ~ moved, dat)
summary(est)
Call:
lm(formula = add_participation ~ moved, data = dat)
Residuals:
Min 1Q Median 3Q Max
-1.788 -1.788 -0.788 0.609 33.609
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.3911 0.1624 8.567 <2e-16 ***
moved 0.3967 0.2010 1.973 0.0488 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.673 on 778 degrees of freedom
(26 observations deleted due to missingness)
Multiple R-squared: 0.00498, Adjusted R-squared: 0.003701
F-statistic: 3.894 on 1 and 778 DF, p-value: 0.04881What if we want to control for the influence of the number of years someone has been at university? Here, we can create a new measure of year and use multiple regression.
#install.packages("modelsummary")
library(modelsummary)
table(dat$class_year)
Year I Year II Year III
320 270 216 dat = dat %>% mutate(
year_bin = as.numeric(case_when(class_year == "Year I" ~ 0,
class_year == "Year II" ~ NA,
class_year == "Year III" ~ 1
)
)
)
table(dat$year_bin)
0 1
320 216 models <- list()
models[['Bivariate']] = lm(add_participation ~ moved, dat)
models[['Multivariate']] = lm(add_participation ~ moved + year_bin, dat)
modelsummary(
models,
estimate = "{estimate}{stars} ({std.error})",
statistic = NULL,
gof_omit = 'IC|RMSE|Log|F|R2$|Std.')| Bivariate | Multivariate | |
|---|---|---|
| (Intercept) | 1.391*** (0.162) | 1.323*** (0.201) |
| moved | 0.397* (0.201) | 0.708** (0.219) |
| year_bin | −0.226 (0.218) | |
| Num.Obs. | 780 | 519 |
| R2 Adj. | 0.004 | 0.021 |
Let’s compare two different ways of defining the year measure.
dat = dat %>% mutate(year_cont = as.numeric(case_when(class_year == "Year I" ~ 0,
class_year == "Year II" ~ 1,
class_year == "Year III" ~ 2 )) )
table(dat$year_cont)
0 1 2
320 270 216 models <- list()
models[['Bivariate']] = lm(add_participation ~ moved + year_bin, dat)
models[['Multivariate']] = lm(add_participation ~ moved + year_cont, dat)
modelsummary(
models,
estimate = "{estimate}{stars} ({std.error})",
statistic = NULL,
gof_omit = 'IC|RMSE|Log|F|R2$|Std.')| Bivariate | Multivariate | |
|---|---|---|
| (Intercept) | 1.323*** (0.201) | 1.537*** (0.202) |
| moved | 0.708** (0.219) | 0.366+ (0.203) |
| year_bin | −0.226 (0.218) | |
| year_cont | −0.146 (0.120) | |
| Num.Obs. | 519 | 780 |
| R2 Adj. | 0.021 | 0.004 |
## Find participation measures that are based on likert
# baseline
bcols = grep("^q13_.*_baseline$", names(dat), value = T)
dat[, paste0(bcols, "_st")] = dat[, bcols]
bcols = paste0(bcols,"_st")
# endline
ecols = grep("^q13_[1-7]_\\d$", names(dat), value = T)
dat[, paste0(ecols, "_st")] = dat[, ecols]
ecols = paste0(ecols,"_st")
# Create treatment variable
dat = dat %>% mutate(moved = case_when(home_region == "Addis Ababa" ~ 0, TRUE ~ 1) )
# clean q13_
levels = c("Never", "Once or Twice", "More than twice", "More than 5 times",
"More than 10 times")
dat = dat %>%
mutate(across(c(bcols),
.fns = ~ factor(.x, levels = levels)))
# Create z-score function from Kling, Liberman, and Katz (2007)
z_score = function(x, y){
# calculate mean and sd of control group
c_mean = mean( as.numeric( unlist(x[, y])) , na.rm = T)
c_sd = sd( as.numeric( unlist(x[, y])) , na.rm = T)
# subtract control group mean; divide by control group SD
( as.numeric(x[, y, drop = TRUE]) - c_mean) / c_sd
}
# calculate z-scores
for (i in c(bcols, ecols)) {
dat[,i] = z_score(dat, i)
}
dat = dat %>%
rowwise() %>%
mutate( z_participation_end = mean(c_across(all_of(bcols)), na.rm = TRUE)) %>%
mutate( z_participation_base = mean(c_across(all_of(ecols)), na.rm = TRUE)) %>%
ungroup()regd = dat %>% select(z_participation_end, z_participation_base, moved, response_id ) %>%
pivot_longer(cols = c(z_participation_end, z_participation_base),
names_to = "time",
values_to = "z_participation") %>%
mutate(time = case_when(time == "z_participation_end" ~ 1,
TRUE ~ 0))
models <- list()
models[['Bivariate']] = lm(z_participation ~ moved, regd)
models[['Multivariate']] = lm(z_participation ~ moved + time, regd)
models[['Interaction']] = lm(z_participation ~ moved + time + moved*time, regd)
modelsummary(
models,
estimate = "{estimate}{stars} ({std.error})",
statistic = NULL,
gof_omit = 'IC|RMSE|Log|F|R2$|Std.')| Bivariate | Multivariate | Interaction | |
|---|---|---|---|
| (Intercept) | −0.159*** (0.026) | −0.156*** (0.030) | −0.156*** (0.037) |
| moved | 0.247*** (0.032) | 0.247*** (0.032) | 0.247*** (0.046) |
| time | −0.004 (0.031) | −0.005 (0.052) | |
| moved × time | 0.001 (0.065) | ||
| Num.Obs. | 1597 | 1597 | 1597 |
| R2 Adj. | 0.035 | 0.034 | 0.034 |